Analysis of B-cell markers in MCS device patients

Statistical Notebook

Introduction

Raw data

The raw data is shown below. There were 196 rows and 42 columns, consisting of 28 patients taken over 7 timepoints. However, there is a significant amount of missing data, resulting in only 120 usable datapoints.

missing.ix <- Reduce(intersect, apply(df[,14:42], 2, function(x) which(is.na(x))))
df.raw <- df[-missing.ix,]
df.raw <- df.raw[order(df.raw$PatientID),]
rownames(df.raw) <- 1:nrow(df.raw)
kable(df.raw, 
      digits = 3,
      row.names = T,
      caption = "Raw Data"
) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%", height = "300px") 

Raw Data

	PatientID	Age	AgeGreater60	Sex	LowIntermacs	InterMACS	RVAD	Sensitized	VAD Indication	Device Type	Survival	Outcome	Time	num Total PBMC	num lymph	lymph	live lymph	CD3 of live lymph	CD19 of live lymph	CD19+CD27-	CD19+CD27+	CD27+38++plasma blasts	CD27-38++ transitional	CD27-IgD+ mature naive	CD27+IgD- switched memory	CD27-IgD- switched memory	CD27+IgD+ unswitched memory	CD27+IgD-IgM+ switched memory	CD27+IgD+IgM+ nonswitched memory	CD19+27+IgG+IgM- memory	CD19+24dim38dim naive mature	CD19+24+38++transitional	CD19CD24hiCD38-memory	CD19+27-38+CD5+transitionals	CD19+CD268+	CD268 of +27-38++transitional	CD19+CD11b+	CD19+CD5+	CD19+CD27+CD24hi	CD19+CD5+CD24hi	CD19+CD5+CD11b+	CD19+27+IgD-38++IgG ASC
1	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	0	169154	35496	20.98	99.53	25.37	19.82	86.02	13.98	1.26	2.54	57.81	11.74	28.06	2.39	21.74	14.36	0.72	81.05	2.86	13.60	0.58	95.22	84.83	10.47	4.33	8.50	1.39	3.13	1.09
2	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	1	63082	9915	15.72	99.26	32.45	26.26	90.06	9.94	0.89	2.09	54.27	7.70	35.71	2.32	31.10	20.47	1.57	84.06	2.59	10.79	0.26	97.02	92.59	7.50	3.13	6.42	1.32	2.44	1.49
3	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	3	75921	21721	28.61	99.31	24.86	33.01	86.38	13.62	1.54	1.26	45.74	10.98	40.44	2.84	31.80	18.04	1.15	86.21	4.82	6.56	0.24	90.41	73.33	10.00	7.25	10.12	3.22	5.67	2.05
4	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	8	213808	36002	16.84	99.19	35.95	14.55	67.94	32.06	4.06	2.17	27.92	28.92	39.64	3.52	20.31	9.10	1.33	68.69	7.10	20.42	1.88	87.38	44.25	10.78	11.12	19.15	5.54	7.35	1.07
5	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	21	106970	31867	29.79	99.18	42.36	15.02	84.69	15.31	1.94	2.11	55.41	13.14	29.13	2.32	20.22	12.88	2.08	83.21	5.20	9.50	1.10	94.44	73.00	9.04	6.99	9.22	2.80	4.78	1.44
6	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	0	119377	77509	64.93	99.61	58.50	7.99	62.67	37.33	9.10	4.82	36.73	31.90	25.62	5.74	11.47	8.95	0.40	69.54	7.69	9.41	2.61	76.44	57.91	23.29	20.76	15.39	5.42	14.90	2.60
7	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	3	96940	57658	59.48	99.70	55.52	9.21	72.15	27.85	3.85	3.82	48.41	21.75	23.49	6.35	12.10	16.04	0.28	75.30	3.72	6.65	2.11	80.52	67.33	15.80	15.23	6.71	1.32	8.86	0.26
8	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	5	154373	73969	47.92	99.65	45.18	8.77	67.36	32.64	6.59	2.77	39.10	26.28	27.95	6.67	17.89	15.71	0.39	71.57	6.26	9.88	1.78	81.86	72.07	20.12	20.17	14.80	5.99	14.65	3.71
9	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	8	155610	63375	40.73	99.51	38.42	12.40	71.69	28.31	4.36	2.62	37.82	23.43	33.50	5.24	17.37	14.92	0.60	65.04	2.80	11.23	2.54	83.55	37.56	15.55	15.45	9.55	2.66	10.91	0.28
10	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	14	244245	115109	47.13	99.47	47.40	9.37	67.50	32.50	3.99	2.08	50.08	26.01	17.22	6.68	15.34	15.25	0.23	77.86	3.21	10.49	1.48	82.13	59.19	21.45	19.77	14.45	2.57	14.94	0.94
11	3	58	younger	Male	Low	3	Yes	No	BTT	HMII	alive	Alive s/p OHT	1	246811	86529	35.06	93.49	67.96	17.48	81.66	18.34	2.05	6.85	51.33	13.69	29.91	5.08	15.68	5.94	0.19	61.95	0.76	13.51	4.07	92.62	34.06	6.90	12.07	8.00	2.54	3.85	0.26
12	3	58	younger	Male	Low	3	Yes	No	BTT	HMII	alive	Alive s/p OHT	3	406771	90343	22.21	97.34	49.28	17.02	75.04	24.96	4.47	18.36	28.45	21.86	46.07	3.62	24.97	4.33	0.13	47.41	1.64	19.70	18.15	72.76	8.08	25.69	28.35	13.95	4.66	19.20	0.41
13	3	58	younger	Male	Low	3	Yes	No	BTT	HMII	alive	Alive s/p OHT	21	330191	81636	24.72	94.57	45.14	14.70	70.51	29.49	8.94	16.23	38.92	23.53	30.83	6.72	23.93	7.07	0.18	46.63	2.35	16.12	8.59	74.45	39.90	24.19	21.18	12.62	3.26	12.93	0.38
14	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	0	309890	41896	13.52	97.03	47.13	8.08	77.31	22.69	2.74	3.29	59.03	19.80	19.37	1.80	21.10	4.46	19.79	71.98	2.41	17.39	0.20	95.80	75.00	5.79	1.89	16.45	0.91	0.67	0.77
15	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	1	274228	15653	5.71	97.85	35.76	28.87	51.45	48.55	2.60	2.13	30.85	43.13	22.77	3.26	20.09	3.92	30.51	57.58	1.31	28.20	0.18	91.23	59.57	12.08	2.04	34.44	0.38	0.16	0.16
16	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	3	591939	57589	9.73	96.76	38.20	42.54	75.71	24.29	2.14	1.33	53.03	21.26	24.16	1.55	25.50	4.02	20.27	67.18	1.14	21.71	0.11	95.28	66.98	3.79	2.43	17.70	0.65	0.84	0.18
17	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	5	189634	46711	24.63	97.74	34.84	37.32	75.30	24.70	0.70	0.55	46.97	22.05	30.28	0.70	23.67	1.56	31.44	73.33	0.75	18.59	0.00	84.61	86.17	3.44	2.11	18.68	0.94	0.63	0.69
18	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	14	182319	40475	22.20	96.28	43.36	11.77	69.34	30.66	1.59	0.89	52.64	27.32	18.56	1.48	22.79	2.88	28.61	75.47	0.33	14.54	0.06	95.90	53.66	7.81	3.27	21.04	1.57	1.50	0.72
19	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	1073919	336538	31.34	100.00	63.29	14.80	93.05	6.95	0.73	8.05	75.49	4.55	19.16	0.80	7.25	8.78	28.04	88.72	2.51	5.18	1.02	90.49	68.18	2.32	4.09	1.69	3.48	1.27	5.78
20	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	1	270301	97229	35.97	91.63	18.43	61.67	97.68	2.32	0.20	5.85	69.45	1.63	28.19	0.73	27.57	21.43	17.05	84.32	7.26	2.73	0.84	99.44	97.60	0.54	2.64	1.51	3.63	0.29	3.40
21	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	3	499016	78019	15.63	89.41	25.83	51.38	92.98	7.02	0.32	4.66	51.34	6.31	41.59	0.76	36.49	5.96	22.81	73.48	6.73	7.93	0.33	98.78	96.52	1.31	2.06	5.19	2.68	0.59	2.45
22	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	115825	32513	28.07	91.04	44.61	15.76	80.52	19.48	4.16	7.82	49.91	18.37	30.52	1.20	14.87	3.33	30.19	61.08	8.70	12.56	1.58	89.80	76.71	2.87	6.07	11.12	8.08	1.59	1.17
23	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	21	407080	97651	23.99	91.29	66.96	12.20	84.66	15.34	2.88	5.60	61.89	13.73	22.65	1.74	20.16	7.23	29.39	68.49	5.28	10.60	1.17	92.50	54.35	3.48	6.53	9.61	8.34	1.85	2.83
24	6	43	younger	Male	High	2	Yes	No	BTT	HMII	alive	Alive s/p OHT	0	515760	266115	51.60	96.30	69.66	8.57	72.22	27.78	6.26	4.77	43.09	26.01	28.26	2.63	16.67	2.67	28.25	50.35	0.41	12.78	4.47	84.76	19.47	9.56	4.55	8.29	0.87	2.18	0.75
25	6	43	younger	Male	High	2	Yes	No	BTT	HMII	alive	Alive s/p OHT	8	539521	223872	41.49	96.72	60.02	12.63	80.13	19.87	1.66	2.48	43.56	19.52	35.24	1.68	20.46	3.67	20.54	55.76	0.18	10.80	2.43	87.82	13.25	6.55	5.96	4.65	0.95	2.51	0.65
26	6	43	younger	Male	High	2	Yes	No	BTT	HMII	alive	Alive s/p OHT	21	564317	216312	38.33	97.40	55.80	4.83	74.19	25.81	5.35	5.77	37.90	25.72	34.76	1.62	17.99	2.18	31.75	61.61	0.77	13.44	6.97	76.55	6.97	10.34	11.14	4.46	0.94	5.14	2.43
27	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	0	89251	65181	73.03	76.76	84.53	1.15	74.65	25.35	1.91	1.39	55.03	22.40	19.62	2.95	15.65	6.80	25.85	70.49	0.87	23.26	0.23	83.16	25.00	8.85	5.03	11.63	1.56	2.95	3.08
28	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	3	375167	234817	62.59	98.46	60.29	10.18	87.72	12.28	0.46	0.48	68.54	9.90	19.03	2.54	30.58	19.69	15.54	72.77	0.13	20.35	0.03	92.56	36.28	3.75	2.70	9.17	1.19	1.27	0.63
29	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	5	256306	162176	63.27	93.82	68.43	8.45	85.45	14.55	0.70	0.53	71.99	11.98	13.37	2.66	32.41	17.96	12.46	69.74	0.41	24.63	0.02	95.72	44.12	4.47	2.23	12.99	1.30	1.20	0.77
30	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	8	463136	330714	71.41	98.79	70.51	7.30	87.51	12.49	1.78	0.58	69.60	10.61	17.82	1.97	29.83	14.55	12.59	75.03	0.15	21.61	0.07	93.59	30.94	4.05	3.07	10.00	1.08	1.35	0.62
31	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	0	92221	27988	30.35	85.45	55.78	1.51	85.56	14.44	0.28	0.28	55.28	4.72	39.72	0.28	0.00	1.89	1.89	89.44	0.00	2.78	0.00	15.00	0.00	3.33	0.83	0.28	0.28	3.33	0.00
32	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	1	273746	46755	17.08	93.25	50.14	7.57	88.27	11.73	0.09	2.12	54.52	3.36	42.00	0.12	8.40	6.87	0.00	92.00	0.12	3.45	0.00	3.94	11.43	1.12	0.64	2.76	0.33	2.15	0.00
33	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	3	156053	71752	45.98	78.38	35.43	27.12	78.57	21.43	0.24	0.60	58.56	5.80	18.33	17.31	12.67	72.96	0.53	94.06	0.49	5.11	0.28	22.79	27.47	0.81	5.53	0.01	0.29	0.19	0.00
34	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	8	156053	71752	45.98	78.38	35.43	27.12	78.57	21.43	0.24	0.60	58.56	5.80	18.33	17.31	12.67	72.96	0.53	94.06	0.49	5.11	0.28	22.79	27.47	0.81	5.53	0.01	0.29	0.19	0.00
35	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	3	116462	41076	35.27	95.44	33.79	25.24	84.51	15.49	2.92	0.99	45.34	13.19	39.12	2.36	17.56	8.75	19.73	79.35	5.17	10.80	0.71	90.04	29.59	10.16	6.81	8.21	3.37	4.86	4.47
36	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	5	222650	85649	38.47	98.40	32.09	19.82	77.54	22.46	5.02	1.10	43.96	19.20	33.50	3.34	17.47	9.77	15.63	71.89	6.68	14.08	0.70	86.49	27.87	13.04	9.84	12.30	4.35	7.98	6.39
37	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	8	102274	63057	61.65	97.73	27.89	20.66	73.75	26.25	7.33	3.45	35.79	21.92	37.73	4.56	9.97	8.20	16.96	63.96	2.40	11.16	3.07	73.97	5.92	16.69	17.77	11.74	4.90	12.19	6.09
38	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	14	133530	81665	61.16	96.46	22.03	16.15	35.30	64.70	13.74	2.33	8.48	52.28	26.37	12.87	5.29	5.79	19.15	47.48	5.65	17.83	5.43	39.88	6.08	52.08	62.11	37.45	21.95	50.15	3.84
39	10	36	younger	Male	High	2	Yes	NA	BTT	PVAD	alive	Alive s/p OHT	0	481928	162052	33.63	99.70	64.73	5.05	3.91	96.09	8.89	1.34	0.72	94.99	3.14	1.15	4.07	0.84	0.32	90.28	2.84	1.58	1.27	6.82	9.17	14.61	90.11	4.95	5.71	13.01	0.51
40	10	36	younger	Male	High	2	Yes	NA	BTT	PVAD	alive	Alive s/p OHT	8	38113	21486	56.37	97.86	58.44	7.88	20.47	79.53	8.09	1.69	1.21	75.00	18.96	4.83	10.07	2.59	0.53	82.19	2.29	9.00	3.59	15.16	0.00	23.07	63.95	4.77	6.52	16.91	0.40
41	10	36	younger	Male	High	2	Yes	NA	BTT	PVAD	alive	Alive s/p OHT	21	80991	19210	23.72	99.00	50.61	3.50	12.33	87.67	14.74	4.81	1.05	85.41	11.13	2.41	8.48	0.52	1.56	80.00	8.27	3.76	7.41	13.68	0.00	30.83	77.74	10.83	14.14	19.85	2.81
42	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	0	411616	119210	28.96	90.82	73.63	5.88	74.30	25.70	2.80	6.46	57.94	19.20	18.72	4.13	19.10	15.42	15.30	53.51	13.78	25.19	0.21	47.23	53.77	5.75	5.92	22.63	3.46	2.55	0.74
43	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	1	430509	79492	18.46	89.31	52.10	20.45	87.13	12.87	1.23	4.76	58.11	9.15	30.55	2.20	29.86	18.05	10.37	63.63	10.23	17.43	0.06	74.05	72.50	2.36	2.98	10.68	1.72	1.03	2.89
44	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	3	395713	119670	30.24	97.90	23.79	6.27	82.17	17.83	0.79	1.99	49.95	13.01	33.91	3.13	28.36	16.83	5.00	56.33	7.13	26.73	0.02	79.50	82.88	1.96	5.55	16.30	3.61	1.88	0.42
45	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	8	329840	73468	22.27	98.63	70.61	3.73	51.50	48.50	8.32	3.77	33.46	41.29	20.92	4.33	25.19	5.89	24.07	46.40	3.73	34.64	0.07	64.14	28.43	10.35	9.61	38.78	4.62	5.55	2.91
46	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	14	350294	88823	25.36	97.63	74.37	4.37	73.35	26.65	4.91	4.83	56.54	21.32	19.84	2.30	34.68	11.25	18.37	61.98	5.78	20.55	0.04	63.91	47.54	6.81	6.17	21.77	2.90	2.74	3.00
47	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	0	128451	60927	47.43	96.82	60.77	4.70	77.26	22.74	2.17	1.95	55.58	18.40	21.54	4.47	27.86	18.89	15.33	69.90	0.94	24.40	0.14	70.73	1.85	8.01	8.66	20.21	5.45	3.32	0.39
48	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	1	142796	22573	15.81	80.67	47.31	2.28	63.37	36.63	9.16	6.51	26.02	33.49	37.11	3.37	38.13	7.50	13.13	53.25	5.78	26.27	1.57	40.24	0.00	27.71	22.17	29.16	18.80	15.66	8.11
49	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	3	157167	38008	24.18	93.47	34.86	3.95	69.71	30.29	3.92	3.92	33.64	26.30	35.78	4.28	39.54	11.49	14.94	66.86	2.85	24.31	1.04	64.50	7.27	21.60	20.53	27.66	14.54	14.18	2.12
50	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	5	121823	23097	18.96	93.48	2.51	2.87	99.52	0.48	0.00	0.00	0.97	0.48	98.55	0.00	50.00	25.00	25.00	94.35	0.00	2.58	0.00	0.16	0.00	17.45	87.40	0.00	0.16	18.26	0.00
51	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	8	196515	88579	45.07	85.54	17.17	3.28	59.38	40.62	3.21	1.41	35.44	35.28	23.87	5.42	62.85	12.32	3.71	66.05	4.38	23.22	0.27	74.53	31.43	25.51	28.12	39.78	25.79	21.57	0.79
52	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	1	87827	15751	17.93	76.38	31.60	12.17	77.46	22.54	6.35	2.53	50.34	20.22	27.05	2.39	19.68	7.62	19.05	75.55	7.38	9.49	0.71	94.26	89.19	11.68	17.96	18.78	6.97	10.86	12.46
53	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	3	55539	13608	24.50	79.00	28.02	15.78	81.90	18.10	6.13	2.42	66.92	14.33	14.98	3.77	20.81	17.45	18.12	78.01	7.08	8.61	0.72	94.04	70.73	11.32	16.69	15.21	5.78	10.50	16.39
54	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	8	295005	36785	12.47	94.58	36.55	9.14	86.32	13.68	2.45	1.82	37.01	12.52	49.18	1.29	16.24	4.71	12.47	8.36	1.73	6.64	0.81	84.78	46.55	12.52	9.28	5.94	1.73	4.47	2.49
55	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	14	60695	20471	33.73	94.91	6.93	7.42	83.97	16.03	3.40	0.69	24.36	14.71	59.33	1.60	18.43	5.07	11.52	54.82	2.91	16.10	0.50	82.44	30.00	51.08	13.25	10.41	3.68	11.31	6.06
56	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	0	52016	30920	59.44	95.77	55.57	24.62	58.73	41.27	2.14	0.97	53.01	5.05	5.45	36.50	9.39	81.29	2.88	81.48	5.65	9.02	0.48	97.75	71.83	0.67	22.89	0.34	1.78	0.49	NA
57	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	1	75516	48801	64.62	92.13	52.90	18.64	70.39	29.61	2.92	2.21	62.18	5.13	7.98	24.70	13.79	69.54	3.19	84.14	3.64	6.36	0.41	94.93	43.78	0.57	15.61	0.23	0.91	0.36	NA
58	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	3	250612	130454	52.05	95.75	60.05	3.29	46.04	53.96	1.34	3.89	33.08	32.47	8.29	26.15	3.56	34.61	0.09	76.71	0.39	13.20	0.28	39.79	14.38	2.26	1.56	10.45	2.89	3.04	0.08
59	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	8	212478	120477	56.70	85.12	48.36	12.51	79.57	20.43	1.35	2.14	58.73	12.87	18.44	9.96	5.63	34.32	0.11	86.74	0.41	6.88	0.15	35.30	53.09	1.23	0.75	2.36	1.04	1.47	NA
60	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	0	86135	70306	81.62	93.71	74.50	1.59	61.74	38.26	1.91	2.67	33.97	14.50	26.81	24.71	20.67	62.47	9.03	73.00	2.00	24.71	1.27	20.23	71.43	0.48	12.79	0.38	3.24	0.48	4.60
61	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	1	81004	41334	51.03	94.06	53.64	2.99	61.70	38.30	3.10	2.75	39.07	13.34	22.12	25.47	12.58	63.11	15.99	78.92	2.93	17.38	3.35	23.24	71.88	0.52	11.70	0.09	1.29	0.52	9.83
62	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	5	90024	39434	43.80	95.32	68.83	5.62	71.59	28.41	0.85	1.66	35.23	13.02	35.80	15.96	22.94	52.37	5.85	76.42	1.14	21.83	0.62	14.54	48.57	0.38	4.69	0.00	0.80	0.38	0.62
63	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	14	106791	80532	75.41	88.82	79.92	1.57	62.32	37.68	1.70	4.02	28.04	22.86	34.11	15.00	23.92	39.41	15.03	61.61	3.66	34.55	1.97	20.36	73.33	0.80	8.57	0.27	1.79	0.98	0.71
64	16	63	older	Male	High	1	No	NA	BTT	HMII	dead	Died	3	366713	77282	21.07	98.61	66.55	12.44	93.65	6.35	0.53	0.60	62.22	4.68	32.38	0.72	19.94	9.89	28.23	83.18	0.51	14.38	0.30	94.72	54.39	2.07	11.33	2.86	3.40	1.48	3.58
65	16	63	older	Male	High	1	No	NA	BTT	HMII	dead	Died	5	262251	55344	21.10	99.55	62.37	14.39	93.00	7.00	1.03	0.69	65.71	4.74	28.37	1.17	25.93	17.28	25.93	79.87	1.20	17.05	0.48	97.76	58.18	2.27	11.13	4.12	4.91	1.87	6.10
66	16	63	older	Male	High	1	No	NA	BTT	HMII	dead	Died	8	309193	61091	19.76	98.15	56.62	7.24	83.61	16.39	3.27	1.52	58.03	11.90	28.13	1.93	24.65	11.36	1.52	79.53	1.45	13.33	0.97	92.63	25.76	3.50	9.23	7.57	3.73	3.48	0.00
67	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	0	147997	54649	36.93	91.75	79.63	7.62	74.86	25.14	0.60	3.45	62.45	22.87	12.30	2.38	45.98	9.28	13.61	74.33	2.67	21.19	0.84	91.31	90.15	3.58	2.25	21.79	1.36	0.58	0.34
68	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	1	177687	56400	31.74	97.13	69.76	12.60	81.83	18.17	0.65	1.61	55.77	16.13	26.02	2.09	53.19	10.70	11.57	82.92	1.38	13.98	0.28	89.83	85.59	3.78	1.85	14.65	1.30	0.68	0.18
69	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	3	233155	67140	28.80	98.44	74.89	5.32	81.96	18.04	1.00	1.88	49.03	16.59	32.90	1.48	49.92	7.20	12.83	83.44	1.39	12.98	0.14	89.70	69.70	4.58	1.62	13.89	1.08	1.08	0.17
70	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	5	99257	29455	29.68	97.33	62.48	4.73	80.60	19.40	1.84	1.25	45.06	18.36	35.25	1.33	44.61	6.32	12.27	82.89	0.88	13.20	0.09	83.55	35.29	7.23	3.54	14.75	2.29	3.02	1.20
71	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	8	76358	18702	24.49	97.32	29.48	11.10	86.29	13.71	2.43	0.50	37.92	12.13	48.32	1.63	46.45	11.70	11.35	83.96	0.20	10.50	0.06	57.82	50.00	9.95	29.60	10.20	0.64	9.36	0.80
72	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	0	56374	12141	21.54	67.38	27.12	3.32	51.84	48.16	1.84	1.84	2.94	45.59	47.79	3.68	0.00	0.00	2.29	92.65	0.37	5.51	0.72	30.15	60.00	6.99	9.56	0.37	0.00	3.68	0.00
73	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	3	7018	3464	49.36	89.03	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
74	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	8	80210	21828	27.21	85.74	20.38	3.50	61.98	38.02	8.55	2.90	5.50	34.96	56.18	3.36	0.00	0.00	3.21	89.31	0.61	2.60	1.49	23.66	0.00	21.37	9.16	1.07	0.15	5.34	1.31
75	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	14	55068	19145	34.77	75.16	16.90	3.30	64.00	36.00	4.21	2.11	10.53	34.11	52.42	2.95	0.00	0.00	2.94	92.42	0.00	3.37	1.00	22.74	0.00	12.42	11.16	0.63	0.00	4.63	0.00
76	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	21	76673	23615	30.80	83.39	26.27	2.06	63.55	36.45	5.42	5.42	5.17	34.48	58.13	2.22	0.00	0.00	2.72	89.41	0.00	3.45	1.95	18.23	4.55	18.23	8.87	1.23	0.25	4.68	0.00
77	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	0	226762	25868	11.41	89.41	25.31	12.05	55.06	44.94	11.31	2.69	14.79	42.57	39.95	2.69	19.87	3.88	32.88	68.95	10.48	13.57	0.59	80.62	22.67	31.91	19.38	15.58	10.37	16.55	8.18
78	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	3	244049	39386	16.14	94.94	40.46	13.23	40.93	59.07	13.36	3.17	8.33	54.90	32.28	4.49	18.52	4.21	23.83	55.18	13.87	19.39	1.44	70.53	35.67	33.62	33.47	22.86	15.26	25.81	5.11
79	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	5	344736	250121	72.55	89.83	7.54	3.81	45.84	54.16	3.96	1.53	13.39	52.70	32.18	1.72	67.14	2.92	2.44	74.68	6.81	13.91	0.21	63.33	58.02	33.05	30.94	19.83	13.96	26.74	0.22
80	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	8	216934	140093	64.58	96.37	4.78	7.22	22.49	77.51	13.49	0.64	7.22	67.61	15.17	10.00	35.88	7.82	17.22	51.28	39.40	6.56	1.33	62.59	58.06	66.67	64.59	53.00	50.26	63.78	4.97
81	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	14	215820	66785	30.94	92.54	17.20	6.01	25.88	74.12	9.15	1.45	3.09	71.72	22.44	2.74	65.43	3.42	0.04	69.38	15.93	8.85	0.74	50.58	24.07	52.78	48.88	32.28	28.19	44.04	0.04
82	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	21	670692	523014	77.98	96.05	5.07	2.64	23.61	76.39	4.11	0.49	4.79	73.53	18.68	3.00	42.40	2.72	4.49	71.95	20.81	5.51	0.45	48.22	44.62	60.13	60.26	39.07	36.82	56.38	1.49
83	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	3	121794	28129	23.10	97.05	52.79	16.42	70.08	29.92	0.29	0.54	50.71	21.46	22.96	4.86	34.18	12.51	23.49	39.78	0.31	55.18	0.35	96.01	16.67	1.74	5.27	23.38	2.95	1.00	0.10
84	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	5	82157	17295	21.05	99.20	72.54	7.41	68.21	31.79	1.10	0.55	48.47	22.90	23.60	5.04	31.90	13.57	18.10	37.84	0.47	56.73	0.35	92.37	14.29	2.05	7.79	26.20	4.72	1.42	1.02
85	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	8	154810	27067	17.48	98.04	72.18	5.03	74.91	25.09	0.90	1.65	58.05	14.76	20.82	6.37	23.43	25.71	22.57	44.94	0.60	49.36	0.71	86.59	22.73	3.30	7.64	20.15	4.87	1.72	0.00
86	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	21	197936	15446	7.80	95.86	50.35	5.44	79.13	20.87	3.98	3.11	52.30	15.90	29.81	1.99	20.23	5.78	24.86	40.62	1.12	50.06	1.89	62.36	8.00	7.20	9.94	13.17	3.60	5.09	1.55
87	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	114586	35996	31.41	100.00	59.03	11.31	75.87	24.13	1.52	7.69	69.83	11.67	10.20	8.30	27.80	29.94	2.04	84.64	0.64	10.54	0.61	91.57	61.34	4.99	7.15	11.57	5.92	2.53	1.92
88	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	3	349054	129419	37.08	95.25	62.35	8.81	79.64	20.36	2.72	4.45	57.70	13.93	25.89	2.48	57.14	13.49	8.69	74.29	3.76	12.35	1.41	85.87	62.53	3.41	5.43	7.67	1.40	1.92	1.06
89	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	5	447309	176789	39.52	93.00	53.57	9.32	80.92	19.08	2.16	5.73	56.72	12.51	28.43	2.34	47.68	12.19	11.00	73.80	4.44	12.53	2.57	83.88	49.09	5.24	6.82	8.32	2.18	2.28	1.72
90	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	558290	301468	54.00	95.62	47.73	7.04	79.13	20.87	3.43	7.49	53.19	13.48	30.32	3.01	36.36	14.54	12.10	70.55	5.15	11.89	3.51	79.79	52.89	6.98	8.50	7.82	1.92	3.09	1.96
91	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	14	479402	362845	75.69	97.95	44.30	6.30	76.53	23.47	3.90	11.10	34.99	12.04	50.31	2.66	21.93	14.88	10.05	70.21	4.66	7.31	7.73	49.96	29.38	13.44	17.37	5.64	3.01	6.95	2.28
92	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	21	334942	191418	57.15	93.28	70.35	4.04	70.89	29.11	6.20	17.70	56.45	19.57	20.44	3.53	38.91	14.64	10.35	55.11	6.61	14.34	8.70	80.96	58.92	7.58	7.96	10.11	2.08	2.08	2.48
93	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	0	129980	41600	32.00	91.40	39.03	2.85	49.63	50.37	7.02	11.55	28.37	29.11	20.79	21.72	19.43	28.34	9.80	69.59	4.90	20.33	6.13	31.15	33.60	20.15	22.27	0.28	4.25	16.54	2.09
94	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	1	59478	20182	33.93	89.07	44.50	10.08	60.93	39.07	2.10	2.70	42.88	13.69	17.66	25.77	13.06	46.12	3.13	81.35	1.16	16.11	0.83	37.58	6.12	5.52	11.64	0.06	1.49	2.76	4.07
95	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	5	26646	9754	36.61	98.69	58.32	7.57	64.88	35.12	2.06	0.55	46.50	12.62	17.70	23.18	15.07	53.31	12.50	83.40	1.10	15.64	1.76	28.12	0.00	7.54	16.46	0.55	1.92	6.04	13.00
96	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	8	73641	21944	29.80	96.79	48.91	4.47	62.74	37.26	2.84	3.16	37.47	17.79	24.95	19.79	19.46	38.11	9.46	79.79	1.89	13.58	1.24	27.05	20.00	14.11	13.37	0.21	2.74	10.53	3.31
97	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	0	44004	20612	46.84	43.65	15.61	21.54	49.64	50.36	10.37	0.57	15.94	45.15	33.18	5.73	24.69	5.68	20.25	44.58	1.55	36.17	0.88	67.23	27.27	32.35	15.48	20.02	4.49	10.94	2.76
98	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	1	99731	23424	23.49	91.81	41.33	10.89	35.24	64.76	11.83	1.62	18.71	59.21	16.19	5.89	28.16	5.89	18.53	39.98	4.66	42.16	2.33	84.07	2.63	22.55	15.21	40.28	6.07	9.70	3.96
99	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	8	129637	50001	38.57	90.65	44.41	6.05	30.66	69.34	17.12	4.05	10.62	66.50	19.60	3.28	26.41	3.46	13.37	39.74	7.15	39.49	11.91	71.97	2.70	32.48	33.14	43.91	7.81	22.70	3.82
100	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	14	134368	46363	34.50	88.32	34.25	4.63	32.47	67.53	12.86	8.38	8.91	64.89	22.88	3.32	24.78	3.89	13.58	36.16	6.01	40.33	13.77	72.11	1.89	29.89	26.36	42.80	6.17	15.97	3.29
101	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	59881	8971	14.98	97.40	50.06	9.57	75.84	24.16	3.59	8.13	64.83	13.40	15.67	6.10	19.31	18.81	16.83	76.67	6.10	6.82	2.00	88.16	47.06	10.53	11.24	10.41	9.93	6.94	11.93
102	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	1	60193	23874	39.66	96.72	58.83	16.24	80.37	19.63	1.87	3.31	68.58	10.62	14.72	6.08	21.33	22.83	3.80	80.61	3.63	8.03	0.33	94.16	57.26	2.96	5.01	10.19	4.72	1.41	1.24
103	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	5	135335	27834	20.57	97.35	44.75	5.61	63.22	36.78	8.49	17.63	39.01	23.16	30.72	7.11	17.53	10.38	12.16	42.70	4.93	19.28	4.68	67.37	4.10	26.64	17.76	19.21	14.54	14.47	10.60
104	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	247113	62042	25.11	96.41	43.93	14.82	78.42	21.58	2.52	5.24	65.35	12.22	17.53	4.90	20.39	17.20	2.93	77.86	2.11	8.78	2.94	84.64	9.70	10.64	13.93	9.90	11.59	8.48	3.51
105	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	14	142053	56983	40.11	93.01	62.72	12.87	85.10	14.90	0.69	1.64	68.02	8.40	21.03	2.55	17.03	13.68	1.08	83.33	0.37	7.38	0.26	79.00	16.07	3.01	5.85	5.16	5.24	1.77	0.70
106	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	21	60384	17029	28.20	93.65	40.96	21.36	83.80	16.20	1.47	2.32	63.52	9.69	23.80	2.99	23.55	15.58	0.91	78.16	1.09	8.13	1.00	83.65	17.72	4.70	8.42	5.99	7.22	3.76	1.56
107	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	47645	37377	78.45	94.61	2.47	9.41	78.56	17.02	1.71	4.63	34.94	16.09	44.77	4.21	4.35	28.99	59.42	67.26	3.10	26.22	7.50	72.82	53.90	13.50	28.20	3.04	7.70	7.67	0.96
108	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	1	87452	43025	49.20	91.54	1.98	26.71	92.06	6.76	0.77	4.67	47.21	4.70	45.17	2.92	6.06	18.18	60.61	81.01	3.95	14.10	12.04	85.97	79.84	5.82	14.15	1.11	4.77	2.84	0.00
109	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	3	124056	80873	65.19	95.41	0.37	14.36	90.47	5.98	0.39	0.60	41.84	4.26	51.08	2.81	18.82	35.29	43.53	86.97	0.10	11.32	1.10	75.49	12.12	3.80	14.94	0.24	0.88	3.25	0.00
110	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	87840	57990	66.02	95.60	0.11	15.80	89.00	9.49	1.56	5.02	45.89	7.18	43.54	3.39	6.10	19.51	56.71	84.23	3.14	9.95	10.70	72.13	52.73	7.06	17.53	0.74	3.13	5.13	0.13
111	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	0	196315	122053	62.17	95.37	60.66	22.66	91.10	8.90	0.45	4.95	85.57	6.30	5.42	2.71	13.13	28.32	20.05	86.32	2.45	5.42	0.96	98.68	94.19	2.50	2.96	4.63	1.24	0.33	0.89
112	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	1	91909	37622	40.93	96.10	43.22	38.85	95.48	4.52	0.14	2.19	85.07	3.09	10.33	1.51	21.03	27.23	9.08	88.84	1.13	3.91	0.23	98.77	88.27	1.16	1.16	2.34	0.38	0.21	0.45
113	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	5	153762	86281	56.11	95.83	65.54	24.24	91.12	8.88	0.65	2.84	81.43	6.30	9.62	2.65	21.14	28.73	15.89	83.00	1.97	6.74	0.36	98.22	93.16	3.20	2.29	5.47	1.11	0.77	1.56
114	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	8	82406	47485	57.62	92.71	65.75	17.46	91.79	8.21	1.08	1.99	76.35	5.70	15.41	2.54	20.56	26.95	11.99	78.67	1.35	7.97	0.28	97.55	85.62	2.65	2.56	5.19	1.18	0.86	1.79
115	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	0	105070	38521	36.66	98.39	47.84	12.88	60.82	39.18	11.08	2.89	33.87	29.43	26.68	10.01	19.72	21.92	2.99	59.35	14.64	20.23	3.45	81.40	46.10	24.47	24.27	24.04	11.00	17.02	1.51
116	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	1	162336	73080	45.02	99.08	49.35	13.95	74.24	25.76	7.77	3.84	45.21	19.33	28.85	6.61	19.16	23.45	0.81	65.34	14.40	13.50	3.17	86.68	71.39	23.26	17.70	15.21	6.93	13.29	0.88
117	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	5	111264	60487	54.36	99.15	45.94	16.13	73.16	26.84	4.64	1.82	52.04	20.07	20.94	6.95	22.93	23.47	0.50	77.11	9.30	9.58	1.20	88.92	81.82	22.20	19.01	16.98	6.89	13.57	1.36
118	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	8	165870	67935	40.96	98.88	43.01	15.86	70.23	29.77	7.85	1.88	40.39	20.49	29.66	9.46	22.70	29.00	1.88	68.61	9.60	16.11	1.47	89.57	69.00	20.12	16.59	17.10	7.02	12.40	2.06
119	28	72	older	Male	Low	4	Yes	NA	BTT	HMII	dead	Died	3	424375	21971	5.18	71.78	56.98	18.53	79.84	20.16	3.08	0.51	57.91	14.07	24.06	3.97	31.74	15.95	13.65	57.94	0.65	28.51	0.17	93.57	53.33	2.87	9.69	15.61	3.87	2.16	3.99
120	28	72	older	Male	Low	4	Yes	NA	BTT	HMII	dead	Died	8	177249	103442	58.36	76.40	15.43	5.64	65.13	34.87	5.21	0.20	35.56	27.91	31.90	4.62	26.60	8.28	13.24	67.89	0.29	20.13	0.00	78.60	22.22	5.50	16.09	25.56	6.80	4.73	0.89

Identification of clusters

We double standardized the data and biclustered it. We found 3 clusters of patients, and 3 clusters of B-cells.

require(pheatmap, quietly = T)
require(RColorBrewer, quietly = T)


colorfun <- function(nlevels, ...){
    if(nlevels > 10) return(standardColors(nlevels))
    if(nlevels <=10) return(pal_d3(...)(nlevels))
} 

annotation.col <- df.raw[,c(1,3,4,6,7,8,9,10,11,12)]
annotation.colors <- lapply(colnames(annotation.col), function(nam){
    colrs <- colorfun(nlevels(annotation.col[[nam]]))
    names(colrs) <- levels(annotation.col[[nam]])
    return(colrs)
})
names(annotation.colors) <- colnames(annotation.col)

df.patient <- double_standardize(t(df.raw[,c(2,14:42)]))
pheatmap(df.patient, 
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
         breaks = c(seq(-max(abs(df.patient), na.rm=T), 0, length.out = 50), seq(0.001, max(abs(df.patient), na.rm=T), length.out = 50)),
         #scale = "row",
         fontsize = 7,
         cutree_rows = 3,
         cutree_cols = 3,
         annotation_col = annotation.col,
         annotation_colors = annotation.colors
)

Metadata associations

We computed the statistical associations between sample groups specified by the metadata.

Survival

We computed metadata factors that were statistically associated with survival.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Survival","Outcome", "LowIntermacs"))], 
                      function(this_factor){
                          table(df[,this_factor], df$Survival)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	5.200	2.650	10.40	1.93e-07	1.11e-06	RVAD
2	NA	NA	NA	3.16e-07	1.11e-06	Device Type
3	NA	NA	NA	0.00153	0.00358	InterMACS
4	2.770	0.931	9.01	0.0513	0.0899	Sensitized
5	1.500	0.797	2.84	0.232	0.273	AgeGreater60
6	0.668	0.328	1.37	0.234	0.273	Sex
7	0.876	0.392	1.89	0.856	0.856	VAD Indication

InterMACS

We computed metadata factors that were statistically associated with a binarized InterMACS score.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("LowIntermacs", "InterMACS"))], 
                      function(this_factor){
                          table(df[,this_factor], df$LowIntermacs)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	2.4e-13	1.92e-12	Device Type
2	NA	NA	NA	2.99e-05	0.000119	Outcome
3	3.190	1.690	6.10	0.00013	0.000346	AgeGreater60
4	1.860	0.965	3.67	0.0493	0.0986	RVAD
5	0.573	0.271	1.21	0.113	0.181	VAD Indication
6	1.330	0.575	3.10	0.557	0.701	Sensitized
7	1.220	0.596	2.47	0.613	0.701	Sex
8	0.955	0.504	1.82	0.88	0.88	Survival

Sensitization

We computed metadata factors that were statistically associated with sensitization.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Sensitized"))], 
                      function(this_factor){
                          table(df[,this_factor], df$Sensitized)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	1.1e-06	5.83e-06	Outcome
2	0.110	0.0347	0.307	1.3e-06	5.83e-06	Sex
3	NA	NA	NA	0.0129	0.0388	InterMACS
4	NA	NA	NA	0.019	0.0427	Device Type
5	0.404	0.1630	0.967	0.0295	0.0531	RVAD
6	2.770	0.9310	9.010	0.0513	0.077	Survival
7	1.330	0.5750	3.100	0.557	0.716	LowIntermacs
8	1.250	0.5310	2.940	0.69	0.776	AgeGreater60
9	0.000	0.0000	Inf	1	1	VAD Indication

Sex

We computed metadata factors that were statistically associated with sex.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Sex"))], 
                      function(this_factor){
                          table(df[,this_factor], df$Sex)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	0.110	0.0347	0.307	1.3e-06	1.17e-05	Sensitized
2	NA	NA	NA	2.89e-06	1.3e-05	Device Type
3	NA	NA	NA	2.21e-05	6.64e-05	Outcome
4	4.030	1.9200	8.870	6.13e-05	0.000138	AgeGreater60
5	NA	NA	NA	0.0309	0.0555	InterMACS
6	1.870	0.7420	5.380	0.227	0.263	VAD Indication
7	0.668	0.3280	1.370	0.234	0.263	RVAD
8	0.668	0.3280	1.370	0.234	0.263	Survival
9	1.220	0.5960	2.470	0.613	0.613	LowIntermacs

Age

We computed metadata factors that were statistically associated with age.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("AgeGreater60"))], 
                      function(this_factor){
                          table(df[,this_factor], df$AgeGreater60)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	9.53e-11	8.58e-10	Outcome
2	NA	NA	NA	5.7e-09	2.56e-08	InterMACS
3	NA	NA	NA	2.35e-08	7.04e-08	Device Type
4	0.216	0.108	0.421	1.25e-06	2.82e-06	RVAD
5	4.030	1.920	8.870	6.13e-05	0.00011	Sex
6	3.190	1.690	6.100	0.00013	0.000195	LowIntermacs
7	1.990	0.931	4.430	0.0577	0.0742	VAD Indication
8	1.500	0.797	2.840	0.232	0.261	Survival
9	1.250	0.531	2.940	0.69	0.69	Sensitized

VAD Indication

We computed metadata factors that were statistically associated with VAD Indication.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("VAD Indication"))], 
                      function(this_factor){
                          table(df[,this_factor], df$`VAD Indication`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	4.18e-30	3.76e-29	Outcome
2	0.000	0.000	0.116	5.89e-10	2.65e-09	RVAD
3	NA	NA	NA	6.72e-05	0.000201	Device Type
4	NA	NA	NA	0.0272	0.0612	InterMACS
5	1.990	0.931	4.430	0.0577	0.104	AgeGreater60
6	0.573	0.271	1.210	0.113	0.17	LowIntermacs
7	1.870	0.742	5.380	0.227	0.292	Sex
8	0.876	0.392	1.890	0.856	0.963	Survival
9	0.000	0.000	Inf	1	1	Sensitized

RVAD

We computed metadata factors that were statistically associated with RVAD.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("RVAD"))], 
                      function(this_factor){
                          table(df[,this_factor], df$RVAD)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	6.86e-20	6.17e-19	Device Type
2	NA	NA	NA	1.39e-13	6.27e-13	Outcome
3	0.000	0.000	0.116	5.89e-10	1.77e-09	VAD Indication
4	5.200	2.650	10.400	1.93e-07	4.35e-07	Survival
5	0.216	0.108	0.421	1.25e-06	2.25e-06	AgeGreater60
6	NA	NA	NA	2.8e-05	4.21e-05	InterMACS
7	0.404	0.163	0.967	0.0295	0.0379	Sensitized
8	1.860	0.965	3.670	0.0493	0.0555	LowIntermacs
9	0.668	0.328	1.370	0.234	0.234	Sex

Device-type

We computed metadata factors that were statistically associated with Device Type.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Device Type"))], 
                      function(this_factor){
                          table(df[,this_factor], df$`Device Type`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = T, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = T, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	1e-04	0.000129	AgeGreater60
2	NA	NA	NA	1e-04	0.000129	Sex
3	NA	NA	NA	1e-04	0.000129	LowIntermacs
4	NA	NA	NA	1e-04	0.000129	InterMACS
5	NA	NA	NA	1e-04	0.000129	RVAD
6	NA	NA	NA	1e-04	0.000129	Survival
7	NA	NA	NA	1e-04	0.000129	Outcome
8	NA	NA	NA	2e-04	0.000225	VAD Indication
9	NA	NA	NA	0.0174	0.0174	Sensitized

Visualization of metadata dependencies

We used the GGally package to visualize pairwise dependencies between sample groups, separating out survivors and non-survivors (indicated by color).

df.temp <- df
colnames(df.temp) <- make.names(colnames(df), unique = T)
#invisible(
suppressWarnings(
    suppressMessages(
        ggpairs(df.temp, 
                mapping = aes(color = Survival), 
                columns = colnames(df.temp)[c(2,4,6,7,8,9,10,11)]) 
    ))

#)

Variability of MCS devices

PCA

Using Principal Component Analysis (PCA), we saw large variability between device types, compared to the variability within device types. We also saw large variability between individual patients. None of the other features were clearly separable.

suppressMessages(require(ggbiplot, quietly = T))
require(ggsci, quietly = T)

isna <- unique(unlist(apply(df[,14:42], 2, function(x) which(is.na(x)))))
pca <- prcomp(double_standardize(df[-isna, 14:42]), center = TRUE, scale. = TRUE)

colorfun <- function(grouping, ...){
    if(nlevels(factor(grouping)) > 10) scale_color_discrete(...)
    if(nlevels(factor(grouping)) <=10) scale_color_d3(...)
} 

plotvars <- names(df)[c(10,1:9,11:17)]
plots.pca <- mclapply(plotvars, function(this_var){
    this_groups <- df[-isna, this_var]
    if(is.numeric(this_groups)) this_groups <- cut(this_groups, 4)
    
    ggbiplot(pca,
             groups = this_groups,
             ellipse = TRUE,
             alpha = 0.3,
             varname.size = 1.2) + 
        colorfun(this_groups, name = this_var) +
        ggtitle(paste0(this_var, " variability")) +
        theme_classic() 
}, mc.cores = detectCores()-1)
names(plots.pca) <- plotvars
#plots.pca$PatientID
#plots.pca$`Device Type`

for(ii in 1:length(plots.pca)){
    cat("  \n###", names(plots.pca)[ii], "\n")
    suppressWarnings(print(plots.pca[[ii]]))
    cat("  \n")
}

Device Type

PatientID

Age

AgeGreater60

Sex

LowIntermacs

InterMACS

RVAD

Sensitized

VAD Indication

Survival

Outcome

Time

num Total PBMC

num lymph

lymph

live lymph

Fisher’s exact test

We analyzed the dependency of device types on the other discrete variables using Fisher’s exact test. We found several dependencies (Benjamini-Hochberg \(q<0.05\)), which we have shown below.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[names(isfactor) != "Device Type"], 
                      function(this_factor){
                          table(df[,this_factor], df$`Device Type`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    fisher.test(this_table, simulate.p.value = T, B = 10000)$p
})
fisher.q <- p.adjust(fisher.p, method = "BH")

signif.ix <- which(fisher.q < 0.05)
signif.order <- sort(fisher.q[signif.ix], index.return = T)$ix
for(this_ix in signif.ix[signif.order]){
    cat("  \n###", names(contingency)[this_ix], "\n")
    cat(paste0("Benjamini-Hochberg qvalue = ", signif(fisher.q[this_ix], 2)),
        ". \n")
    print(kable(contingency[this_ix][[1]],  row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
}

AgeGreater60

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
younger	56	0	7	21	7
older	84	7	7	0	7

Sex

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
Female	35	0	0	14	0
Male	105	7	14	7	14

LowIntermacs

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
Low	63	0	14	0	0
High	77	7	0	21	14

InterMACS

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
1	7	0	0	0	7
2	70	7	0	21	7
3	56	0	14	0	0
4	7	0	0	0	0

RVAD

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
No	105	7	14	0	0
Yes	35	0	0	21	14

Survival

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
alive	98	0	14	7	7
dead	42	7	0	14	7

Outcome

Benjamini-Hochberg qvalue = 0.00013 .

	HMII	CMAG	HVAD	PVAD	TAH
Alive	28	0	0	0	0
Alive s/p OHT	70	0	14	7	7
Died	35	7	0	7	0
Died post OHT	7	0	0	0	0
Died s/p OHT	0	0	0	7	7

VAD Indication

Benjamini-Hochberg qvalue = 0.00045 .

	HMII	CMAG	HVAD	PVAD	TAH
BTT	98	7	14	21	14
DT	42	0	0	0	0

Sensitized

Benjamini-Hochberg qvalue = 0.018 .

	HMII	CMAG	HVAD	PVAD	TAH
No	42	0	7	0	7
Yes	28	0	7	7	7

One-way repeated measures anova

We attempted to analyze the differences in B-cell levels across device types using a linear mixed effect model. Here we report variables that had a statistically significant variance (Benjamini-Hochberg \(q<0.05\)) across devices, time, or their interaction.

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- colnames(df)[14:42]
colnames(pvals) <- c("pvalue")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- c("qvalue")

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)

sig.qtable <- cbind(pvals,qBH)[sigvars,,drop=F][order(apply(pvals[sigvars,,drop=F], 1, min)),,drop=F]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")

Significant device-type q-values

	pvalue	qvalue	B-cell
1	0.011	0.31	CD27-IgD+ mature naive
2	0.045	0.45	CD19+CD5+

We plotted the average across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the fit. Note that the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- sumtable$coefficients[-1, c(1:5)] # drop intercept
        sigsum <- sumtable[sumtable[,5] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,5]),,drop=F], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD27-IgD+ mature naive

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD	-33.85486	10.06687	23.54824	-3.362998	0.0026287

CD19+CD5+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD	23.16366	8.615765	23.18582	2.688521	0.0130627

Mixed-effect model

We attempted to analyze the differences in B-cell levels across device types using a linear mixed effect model. We used time as a continuous variable, and included the interaction term. Here we report variables that had a statistically significant variance (Benjamini-Hochberg \(q<0.05\)) across devices, time, or their interaction.

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type * Time + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1,2,3)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- varnames
colnames(pvals) <- c("Device", "Time", "Device:Time")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- colnames(pvals)

sigvars <- apply(apply(qBH, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)
sig.qtable <- qBH[sigvars,][order(apply(qBH[sigvars,], 1, min)),]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        Device = cell_spec(Device, color = ifelse(qtable$Device > 0.05, "grey", "red")),
        `Device:Time` = cell_spec(`Device:Time`, color = ifelse(qtable$`Device:Time` > 0.05, "grey", "red")),
        `Time` = cell_spec(`Time`, color = ifelse(qtable$`Time` > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")

Significant device-type q-values

	Device	Time	Device:Time	B-cell
1	0.6	0.0023	0.007	CD19+CD11b+
2	0.5	0.025	0.82	CD268 of +27-38++transitional
3	0.82	0.025	0.82	CD19+CD5+CD11b+

We plotted the average across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit. Note that the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- sumtable$coefficients[-1, c(1:5)] # drop intercept
        sigsum <- sumtable[sumtable[,5] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,5]),,drop=F], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD19+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:Time	2.1101624	0.5538816	89.91928	3.809771	0.0002541
Device.TypePVAD:Time	0.8054722	0.2943289	88.29844	2.736640	0.0075041
Time	0.2628129	0.1188074	89.80591	2.212092	0.0294973

CD268 of +27-38++transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	-0.7873216	0.3030382	89.7994	-2.598094	0.0109556

CD19+CD5+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.290629	0.1150892	91.27507	2.525249	0.0132869

Two-way repeated measures anova

We attempted to analyze the differences in B-cell levels across device types using a mixed effect model (a.k.a. two-way repeated measures ANOVA). Here we report variables that had a statistically significant variance (Benjamini-Hochberg \(q<0.05\)) across devices, times, or their interaction. As there were 5 device types and 7 timepoints, but only 120 samples, this model is severely underpowered. The posterior belief in any of these results should therefore be quite small (as a consequence of Bayes rule).

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type * factor(Time) + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1,2,3)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- varnames
colnames(pvals) <- c("Device", "factor(Time)", "Device:factor(Time)")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- colnames(pvals)

sigvars <- apply(apply(qBH, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)
sig.qtable <- qBH[sigvars,][order(apply(qBH[sigvars,], 1, min)),]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        Device = cell_spec(Device, color = ifelse(qtable$Device > 0.05, "grey", "red")),
        `Device:factor(Time)` = cell_spec(`Device:factor(Time)`, color = ifelse(qtable$`Device:factor(Time)` > 0.05, "grey", "red")),
        `factor(Time)` = cell_spec(`factor(Time)`, color = ifelse(qtable$`factor(Time)` > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")

Significant device-type q-values

	Device	factor(Time)	Device:factor(Time)	B-cell
1	0.16	0.071	3e-10	CD19+CD5+
2	0.33	0.041	8.7e-06	CD27+IgD- switched memory
3	0.54	0.071	7.8e-05	CD19+CD27+
4	0.54	0.31	7.8e-05	CD27-IgD- switched memory
5	0.59	4e-04	0.0017	CD19+CD11b+
6	0.071	0.16	0.0034	CD27-IgD+ mature naive
7	0.93	0.0058	0.98	CD19 of live lymph
8	0.93	0.94	0.012	CD19CD24hiCD38-memory
9	0.27	0.041	1	CD27+IgD-IgM+ switched memory

We plotted the average across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit. Note that the reference level for the time comparisons is timepoint 0, and the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- sumtable$coefficients[-1, c(4:5)] # drop intercept
        sigsum <- sumtable[sumtable[,2] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,2]),], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD19+CD5+

	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	7.156565	0.0000000
Device.TypePVAD	3.590703	0.0008326
Device.TypePVAD:factor(Time)3	-2.877181	0.0052873
Device.TypePVAD:factor(Time)5	-2.809579	0.0063611
Device.TypePVAD:factor(Time)8	-2.720946	0.0082277

CD27+IgD- switched memory

	t value	Pr(>|t|)
Device.TypePVAD:factor(Time)3	-5.406064	0.0000008
Device.TypePVAD	4.233715	0.0001521
Device.TypePVAD:factor(Time)8	-3.571943	0.0006518
Device.TypePVAD:factor(Time)5	-3.404439	0.0010914
Device.TypeHVAD:factor(Time)1	2.457443	0.0165160
Device.TypeTAH:factor(Time)1	2.120071	0.0376068

CD19+CD27+

	t value	Pr(>|t|)
Device.TypePVAD:factor(Time)3	-4.796905	0.0000087
Device.TypePVAD	3.502922	0.0011749
Device.TypePVAD:factor(Time)5	-3.024465	0.0034502
Device.TypePVAD:factor(Time)8	-2.789363	0.0068193
Device.TypeHVAD:factor(Time)1	2.452316	0.0167262

CD27-IgD- switched memory

	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	5.189516	0.0000021
Device.TypeTAH:factor(Time)14	2.262837	0.0264767
Device.TypePVAD:factor(Time)3	-2.144030	0.0352412

CD19+CD11b+

	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)14	5.250946	0.0000015
Device.TypeTAH:factor(Time)1	2.722190	0.0081768
Device.TypePVAD:factor(Time)14	2.473834	0.0157228
Device.TypeTAH:factor(Time)3	2.083219	0.0408780

CD27-IgD+ mature naive

	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	-3.905748	0.0002179
Device.TypePVAD	-3.456469	0.0012664
Device.TypeTAH:factor(Time)14	-2.953963	0.0042493
Device.TypeTAH:factor(Time)1	-2.579593	0.0120122
factor(Time)8	-2.477197	0.0157031
Device.TypeHVAD:factor(Time)1	-2.083232	0.0409155

CD19 of live lymph

	t value	Pr(>|t|)
factor(Time)1	3.054598	0.0031913
Device.TypeHVAD:factor(Time)5	2.297061	0.0246707
Device.TypeHVAD:factor(Time)3	2.039507	0.0452529

CD19CD24hiCD38-memory

	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	-3.776490	0.0003363
Device.TypeHVAD:factor(Time)1	2.344675	0.0219407
factor(Time)1	-2.161244	0.0341830

CD27+IgD-IgM+ switched memory

	t value	Pr(>|t|)
factor(Time)5	3.237317	0.0018346
factor(Time)3	2.819578	0.0062189
factor(Time)14	2.232742	0.0287354

HeartMate-II analysis

The HeartMate-II (HMII) recipients were the largest group, and we analyzed them by themselves due to the previously observed variability across devices.

require(reshape2, quietly = T)
df.HMII <- subset(df, df$`Device Type`=="HMII")
df.long <- melt(df.HMII, id.vars = colnames(df)[1:13])

groups <- make.names(c("AgeGreater60", 
                       "Sex",
                       "LowIntermacs",
                       "RVAD", 
                       "Sensitized",
                       "VAD Indication", 
                       #"Device Type", 
                       "Survival",
                       "Outcome"))

names(df.long) <- make.names(names(df.long))

plots.ts <- mclapply(unique(df.long$variable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            scale_color_aaas() + scale_fill_aaas() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- unique(df.long$variable)

# for(ii in c(1:length(plots.ts))){
#     for(jj in 1:length(plots.ts[[ii]])){
#         suppressWarnings(print(plots.ts[[ii]][[jj]]))
#     }
# } 

PCA

Using PCA, we found large variability between individual patients, compared to the variability within individual patients. None of the other features were clearly separable.

suppressMessages(require(ggbiplot, quietly = T))
require(ggsci, quietly = T)

# Efron's double standardization
double_standardize <- function(x, niter = 100) {
    for(i in 1:niter) x <- t(scale(t(scale(x))))
    return(as.data.frame(x))
}

isna <- unique(unlist(apply(df.HMII[,14:42], 2, function(x) which(is.na(x)))))
pca <- prcomp(double_standardize(df.HMII[-isna, 14:42]), center = TRUE, scale. = TRUE)

colorfun <- function(grouping, ...){
    if(nlevels(factor(grouping)) > 10) scale_color_discrete(...)
    if(nlevels(factor(grouping)) <=10) scale_color_d3(...)
} 

plots.pca <- mclapply(names(df.HMII)[1:17], function(this_var){
    this_groups <- df.HMII[-isna, this_var]
    if(is.numeric(this_groups)) this_groups <- cut(this_groups, 4)
    
    ggbiplot(pca,
             groups = this_groups,
             ellipse = TRUE,
             alpha = 0.3,
             varname.size = 1.2) + 
        colorfun(this_groups, name = this_var) +
        ggtitle(paste0(this_var, " variability")) +
        theme_classic() 
}, mc.cores = detectCores()-1)

names(plots.pca) <- names(df.HMII)[1:17]
#plots.pca$PatientID
#plots.pca$`Device Type`

for(ii in 1:length(plots.pca)){
    cat("  \n###", names(plots.pca)[ii], "\n")
    suppressWarnings(print(plots.pca[[ii]]))
    cat("  \n")
}

PatientID

Age

AgeGreater60

Sex

LowIntermacs

InterMACS

RVAD

Sensitized

VAD Indication

Device Type

Survival

Outcome

Time

num Total PBMC

num lymph

lymph

live lymph

Two-way repeated measures anova

We attempted to analyze the differences in B-cell levels for various features using a mixed effect model (a.k.a. two-way repeated measures ANOVA). Here we report variables that had a statistically significant variance (Benjamini-Hochberg \(q<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " * factor(Time) + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1,3)]))


# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.ranked <- qvals.df[order(qvals.df$qvalue, decreasing = F),]
shortlist <- qvals.df.ranked[qvals.df.ranked$qvalue <= 0.3,]


kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")

Significant results

	B-cell	parameter	qvalue
2	num lymph	AgeGreater60	0.040
145	CD19+27+IgD-38++IgG ASC	LowIntermacs	0.040
52	CD268 of +27-38++transitional	AgeGreater60:factor(Time)	0.073
276	CD27+IgD-IgM+ switched memory	Sensitized:factor(Time)	0.081
61	lymph	Sex	0.200
348	CD19+27+IgD-38++IgG ASC	VAD.Indication:factor(Time)	0.213
232	CD19+27+IgD-38++IgG ASC	RVAD:factor(Time)	0.230
155	CD27-38++ transitional	LowIntermacs:factor(Time)	0.254
3	lymph	AgeGreater60	0.257

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- sumtable$coefficients[-1, c(1:5)] # drop intercept
    sigsum <- sumtable[sumtable[,5] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,5]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

num lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)21	169798.9	44963.12	1780.755	3.776403	0.0001643
AgeGreater60older:factor(Time)21	-214015.6	62340.16	1651.364	-3.433029	0.0006116
AgeGreater60older:factor(Time)5	-136885.7	61168.57	2679.738	-2.237844	0.0253131
factor(Time)14	128214.1	57636.41	1470.433	2.224532	0.0262639

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
LowIntermacsHigh	3.095	1.279721	69	2.418497	0.0182277

CD268 of +27-38++transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
AgeGreater60older:factor(Time)8	-27.22177	10.33554	54.58048	-2.633803	0.0109633
factor(Time)14	-24.92941	12.43679	55.44507	-2.004489	0.0499111

CD27+IgD-IgM+ switched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SensitizedYes:factor(Time)3	17.965877	5.302290	23.31974	3.388324	0.0024949
factor(Time)5	9.747852	4.547781	23.30369	2.143431	0.0427391

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-22.22807	9.001586	53.16972	-2.46935	0.0167811

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)5	4.722429	1.330913	60.14685	3.548263	0.0007592
VAD.IndicationDT:factor(Time)5	-6.577445	1.874482	57.09611	-3.508940	0.0008850

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
RVADYes:factor(Time)5	11.55470	2.912778	59.01744	3.966900	0.0001997
factor(Time)14	-2.17979	1.043488	56.60858	-2.088945	0.0412202

CD27-38++ transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
LowIntermacsHigh:factor(Time)21	-11.380113	2.894440	59.65377	-3.931715	0.0002222
factor(Time)21	9.440964	2.450218	61.33982	3.853112	0.0002815

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)